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ABSTRACT 

Numerous formulations of finite volume schemes for the Euler and Navier- Stokes equations 
exist, but in the majority of cases they have been developed for structured and stationary 
meshes. In many applications, more flexible mesh geometries that can dynamically adjust to 
the problem at hand and move with the flow in a (quasi) Lagrangian fashion would, how- 
ever, be highly desirable, as this can allow a significant reduction of advection errors and an 
accurate realization of curved and moving boundary conditions. Here we describe a novel 
formulation of viscous continuum hydrodynamics that solves the equations of motion on a 
Voronoi mesh created by a set of mesh-generating points. The points can move in an arbitrary 
manner, but the most natural motion is that given by the fluid velocity itself, such that the 
mesh dynamically adjusts to the flow. Owing to the mathematical properties of the Voronoi 
tessellation, pathological mesh-twisting effects are avoided. Our implementation considers 
the full Navier-Stokes equations and has been realized in the AREPO code both in 2D and 3D. 
We propose a new approach to compute accurate viscous fluxes for a dynamic Voronoi mesh, 
and use this to formulate a finite volume solver of the Navier-Stokes equations. Through a 
number of test problems, including circular Couette flow and flow past a cylindrical obstacle, 
we show that our new scheme combines good accuracy with geometric flexibility, and hence 
promises to be competitive with other highly refined Eulerian methods. This will in particu- 
lar allow astrophysical applications of the AREPO code where physical viscosity is important, 
such as in the hot plasma in galaxy clusters, or for viscous accretion disk models. 

Key words: hydrodynamics - methods: numerical. 
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1 INTRODUCTION 

The last two decades have seen remarkable advances in the nu- 
merical solution of the compressible Navier-Stokes (NS) equations, 
which lies at the heart of computational fluid dynamics (CFD) and 
computational aeroacoustics, but also as numerous applications in 
astrophysics. In particular, important progress has been made in 
approaches based on the finite volume method (FVM), both using 
structured as well as unstructured grids (see Mavriplis 1997, for 
a review). Other popular techniques include finite element meth- 
ods (FEM), discontinuous Galerkin schemes, and even mesh-free 
approaches such as smoothed particle hydrodynamics (Sijacki & 
|Springel|2006l ). 

When unstructured grids have been employed, they were most 
most often in the form of triangular grids in two dimensions, or 
tetrahedral grids in three dimensions. Indeed, finite- volume im- 
plementations of the two-dimensional NS equations on triangular 
meshes date back to work by Mavriplis & Jameson (1990), Frink 
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(1994} and |CoirieF& Powell (1996). Much recent work has also 
focused on developing optimum mesh-generating algorithms that 
require minimal human input and yield efficient representations of 
geometrically complex simulation domains. However, little work 
has been done on dynamically evolving meshes, such as those we 
shall consider here. 

Because unstructured meshes have been demonstrated to be 
accurate and efficient for both steady-state and transient compress- 
ible inviscid flows (Barth 1992, Venkatakrishnan 1996), they are 
now used regularly in engineering applications. Moreover, the ge- 
ometric flexibility of unstructured grids allows the use of simple 
coordinate systems (in the laboratory frame) without the need to 
work with complex coordinate transformations to describe curved 
surfaces (e.g. see Toro 2009 ). Indeed, hard boundaries can be tai- 
lored by carefully positioning a few cell faces or mesh generating 
points along the surface, and creating the triangulation through De- 
launay tessellation. As a result, most NS applications on unstruc- 
tured meshes for industrial design make use of triangular grids, typ- 
ically based on the finite element method, although finite volume 
schemes have also been considered. Detailed reviews and stability 
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analysis of explicit FVM for the NS equations on Cartesian and De- 
launay meshes can be found, e.g, in the doctoral theses of Coirier 
( fT994] > and |Munito-ishna| ( [2009l >. 

In this work, we present a numerical scheme that solves the 
NS equations on a general unstructured moving mesh that is con- 
structed as the Voronoi tessellation of a distributed set of points 
that move with the local velocity field. Despite being, in the general 
sense, an "unstructured" mesh, the Voronoi diagram has a mathe- 
matically well-defined structure that makes the resulting schemes 
comparatively simple and robust (e.g.[Mishev 1998). In fact, this 
type of mesh is commonly adapted for the construction of finite 
volume methods for elliptic problems and has been in use in nu- 
merical studies of solid state physics (Sukumar et al. 1998, Suku- 
mar 2009) such as simulations of fractures and cracks (Sukumar & 
Bolander 2009 ), as well as numerical simulations of oil reservoirs. 
Some studies (Christov 2009) have also examined how reconstruc- 
tions designed for unstructured triangulations can be extended to 
static Voronoi meshes. 

However, Voronoi meshes have infrequently been applied to 
hyperbolic conservation laws such as the Euler equations, let alone 
moving Voronoi meshes. To our knowledge, the earliest attempts 
to use dynamically adaptive Voronoi tessellations for the NS equa- 
tions data back to Borgers & Peskin] ([1987 ), although for very 
simplified, incompressible, two-dimensional problems. Around the 
same time, Dukowi cz et al.| (T989) developed the General Topol- 
ogy Godunov Method. This method - based on a mesh that is not 
quite a Voronoi tessellation, but similar in spirit - was introduced as 
an alternative to the Lagrangian particle methods (see, for example 
Brackbill & Monaghan 1988 ) which gained increasing popularity 
in computational plasma physics and astrophysics in the following 
years. 

Recently, a complete three-dimensional implementation of the 
Euler equations on a moving Voronoi mesh has been described and 
implemented in the AREPO code by Springel (2010) (see also Duf- 
fell & MacFadyen 201 1). The work we present here is an extension 
of the AREPO scheme to the NS equations, which we have real- 
ized in this code as an optional module. AREPO can be classified 
as an arbitrary Lagrangian/Eulerian (ALE; |Hirt et al.|1974| > code, 
in the sense that the mesh can be moved with the velocity of the 
flow so that quasi-Lagrangian behavior results and the mass flux 
between cells is minimized (although it is not strictly zero, in gen- 
eral). On the other hand, the mesh may also be kept stationary if 
desired, effectively yielding an Eulerian formulation. We note that 
because the mesh-generating points may also be arranged on a reg- 
ular lattice and arbitrarily refined with time, the AREPO code natu- 
rally includes ordinary Eulerian techniques on a Cartesian grid and 
adaptive mesh refinement (AMR) algorithms as special cases. 

Besides the work of Duffell & MacFadyen (2011), the new 
Voronoi- ALE method of Norris et al. (2010), which includes vis- 
cous terms, is the approach most closely related to that presented 
here, although it is restricted to the incompressible NS equations. 
Also, |Ata et al.| ( [2009] > have applied a Voronoi-based finite volume 
scheme to the two-dimensional inviscid shallow water equations, 
in terms of an algorithm they referred to as the 'natural volume' 
method. 

Although primarily designed for astrophysical fluid dynam- 
ics where self-gravity is an important ingredient (see for example 
Vogelsberger et al. 2011), the moving Voronoi mesh approach of 
AREPO offers a number of features than can be advantageous for 
more general problems in fluid dynamics. First, the moving mesh 
geometry is adaptive in a continuous manner and can naturally re- 
spond to the local flow, increasing the resolution automatically and 



smoothly in regions where the flow converges. (In contrast, AMR 
codes refine the grid discontinuously in time, which can introduce 
errors that are potentially difficult to assess.) Importantly, this La- 
grangian character of the dynamics yields reduced advection er- 
rors and a very low numerical diffusivity of the scheme. Second, 
the moving mesh formulation retains the Galilean-invariance of the 
fluid dynamics at the discretized level of the equations (Springel 
|2010| ). In other words, the truncation error of the scheme does not 
depend on the bulk velocity of the system, unlike for traditional 
Eulerian and AMR codes, and the quality of the solution does not 
degrade when high-speed flows are present. While conventional 
fixed-mesh Eulerian codes may, in principle, be able to suppress 
additional errors from large bulk velocities by using a sufficiently 
fine mesh (see Robertson et al. 2010 for a study of Galilean in- 
variance in grid codes), this strategy can become computationally 
prohibitive, and it also depends on the magnitude of the bulk veloc- 
ity involved. It is therefore desirable to construct efficient methods 
that yield manifestly Galilean-invariant solutions (modulo floating 
point round-off errors). Third, the moving mesh approach allows 
much larger timesteps in the case of rapidly moving flows, because 
it can avoid the At < d/v stability constraint (where d is the cell 
size and v the bulk velocity) that augments the Courant condition 
in the Eulerian case. 

From an astrophysical standpoint, compressible viscous flow 
remains a viable approximation to more complex or computation- 
ally expensive momentum transport mechanisms such as magneto- 
hydrodynamic turbulence or anisotropic plasma viscosity. Global 
simulations of cold accretion disks around protostellar objects (e.g. 
see |de Val-Borro et al. 2006 ) still include shear viscosity coeffi- 
cients in the form of a Shakura-Sunyaev eddy viscosity coefficient 
( |Shakura & Sunyaev|1973] >. 

An even clearer case for the need of a viscous treatement of 
astrophysical gasdynamics is given by the interacluster medium of 
hot galaxy clusters. Here the Spitzer-Braginskii viscosity (Bragin- 
skii 1965 ) becomes quite significant, certainly in the unmagnetized 
case, which has been studied both using grid (Ruszkowski et al. 
|2004] ) and SPH jSijacki & Springel||2006] ) codes. In this regime, 
the commonly adopted assumption of inviscid behaviour with an 
effectively infinite Reynolds number is in principle incorrect and 
should in future simulation work be replaced with a full accounting 
of the correct physical viscosity. 

Additionally, physical viscosity can be implemented on tur- 
bulent cascades with resolved inertial range (see Bauer & Springel 
|2011| for an application of our viscosity approach) in order to pre- 
scribe a well-specified Reynolds number and a physically correct 
shape for the dissipation range, unaffected by the details of the 
numerical viscosity of the hydro scheme, which would otherwise 
induce the dissipation of turbulence on the grid scale. This can in 
particular inform the ongoing debate whether artificial viscosity ef- 
fects in SPH can affect the turbulent cascade above the formal res- 
olution length ( Bauer & Springel|2011[|Price|2012| ). 

This paper is organized as follows. In Section 2, we briefly 
review the basic NS equations we want to solve, and the role 
and meaning of the different viscosity coefficients. In Section 3, 
we then introduce in detail our discretization and time integration 
schemes, emphasizing a description of the calculation of suitable 
velocity gradient estimates at face centers, and of second-order 
derivatives of the velocity field. We then move on to discuss the 
performance of our new approach for a number of test problems in 
Section 3. Finally, we summarize our results and present our con- 
clusions in Section 4. 
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2 THE NAVIER-STOKES EQUATIONS 

The compact form of the Euler equations, when written in terms of 
the vector of conserved quantities U (Toro 2009 ) is 

d t U + V • F adv (U) - 0, 



with 
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is the mass-momentum-energy flux density tensor (3 x 5). The op- 
erator V • ( ) in Eq. |TJ is a tensor divergence, i.e. in tensor nota- 
tion we have {V • F a d v } a = #6 -Fadv ha '. The momentum compo- 
nents in the conservative form of Equation |TJ represent a transfer 
of momentum, owing merely to the mechanical transport of differ- 
ent particles of fluid from place to place and to the pressure forces 
acting on the fluid (e.g. Landau & Lifshitz 1959). In Eq. {I} we 
have made the advective character of the fluxes explicit by denot- 
ing them Fadv 

The internal friction present in any real fluid causes an irre- 
versible transfer of momentum from points where the velocity is 
large to those where it is small. The momentum flux density tensor 
is thus altered from its ideal from in Eq. |5]), where it only contains 
an inertial and an isotropic component (described by a symmetric 
stress tensor due to the local pressure P), to a modified expression 
that accounts for an irreversible viscous transfer of momentum 



v + pi - n, 



(4) 



pv v + fl — > pv 

where PI — II is the total stress tensor and II is called the viscous 
stress tensor. The latter includes the effects of isotropic compres- 
sion and expansion forces ("bulk viscosity") as well as shearing 
forces ("shear viscosity"). 

Similarly, the energy component of Eq. (3} is affected by the 
inclusion of the viscous stress tensor. Because of the dissipative 
nature of viscosity, a conservative formulation of the NS equations 
must include a contribution of II to the energy budget, i.e. the work 
per unit area per unit time, 



(pe + P) v — ► (pe + P) v - IIv 



(5) 



needs to explicitly account for the work done by viscous forces. 

A general parametrization of the viscous stress tensor II is 
given by 



11 = 77 



[Vv + (Vv) T ] 
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v) +CI(V-v) 



(6) 



Often, the viscous stress tensor is decomposed into a traceless part 
and a diagonal part, such that the first corresponds to constant- 
volume shear deformations (often called the rate -of- deformation 
tensor) and the second to isotropic expansions/contractions. Ac- 
cordingly, r\ in Eq. {6} is commonly referred to as the shear vis- 
cosity and C as the bulk viscosity. The degree of resistance to uni- 
form contractions/expansions is intrinsic to the molecular/chemical 



properties of the fluid in question, and can be understood through 
kinetic theory. In this picture, bulk viscosity arises because kinetic 
energy of molecules is transferred to internal degrees of freedom. 
Ideal monoatomic gases (modeled as hard spheres interacting only 
through elastic collisions) have no internal degrees of freedom, and 
are thus expected to have vanishing bulk viscosity. At one time 
Stokes suggested that this might in general be true (the so-called 
Stokes ' hypothesis of £ = 0) but later wrote that he never put much 
faith in this relationship (Graebel 2007). Indeed, when deviations 
from the ideal gas equation of state are included in a hard- sphere, 
Chapman-Enskog approach to kinetic theory, a non-zero value for 
the bulk viscosity is obtained. In an extension of the hard sphere 
fluid model, the Longuet-Higgins-Pople relation ( — (5/3)7? re- 
sults (March 2002 ), motivating the hypothesis that both viscosities 
are always related in a linear fashion (but see Meier et al. 2005 ). 
In general, we consider r\ and ( as essentially arbitrary input prop- 
erties to our simulations, which may also depend on local physical 
parameters such as temperature or density. Although the effects of 
physical bulk viscosity are not harder to implement numerically 
than those of shear viscosity, the physical origin of bulk viscos- 
ity is often less clear. Also, we note that many numerical solvers 
for viscous flow focus on the incompressible regime (V • v = 0), 
where the existence of a physical bulk viscosity is in any case not 
of importance. However, for compressible flow, the value of £ may 
still become important in certain situations. 

When the effects of viscosity are included, the formerly ho- 
mogeneous differential equations of the Euler form (Eq.[TJ become 



! tU + V-F adv (U) = S(U) 



where S (U) is a viscous source term given by 



(7) 



S(U) 



o , n , nv 



(8) 



The solution of the Euler equations with source terms is often 
handled by operator- splitting methods (e.g. [Toro 2009, LeVeque 
2002 ). That is, the numerical scheme alternates between an advec- 
tive step that solves the homogeneous part, and a source- term step. 
Thus, the solution of Eq. ([7]) is split into a two stage problem: 
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(9) 
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Typically, the source terms are more easily written in the prim- 
itive variable formulation of the Euler equations. A common 
choice of the primitive-variable vector is W = (p,v,P) T — 
(p, v x , v y , v z , P) T , which we also adopt here. For sources corre- 
sponding to the NS viscous terms (Eq.|S]), only the v component of 
W is affected, thus simplifying the solution method of the source- 
term step. The three-dimensional Euler equations can be written in 
the primitive variable form as (Toro 2009) 



d t W + Ai(W) d x W + A 2 (W) d y W + A 3 (W) d z W = 0. 

(11) 
For this choice of variables, the coefficient matrices are given by 
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which is exactly equivalent to the familiar equations 
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dt ' '~ dxi ' " l dx, 
In this formulation, the viscous terms of the NS equations, 
affect only the velocity, are (e.g. |Landau & Lifshitz 1959| > 

S(W) = if r,Av + (C + §»?) V (V • v) J . 



(15a) 

(15b) 

(15c) 
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(16) 



An alternative to expressing the viscosity effects as source 
terms is to absorb them directly into the flux divergence, 



'*U + V 



(U) 



(U) 



= 0, 



(17) 



which highlights the still conservative character of the NS equa- 
tions. Here diffusive fluxes, defined by 



Fciff (u) = o , n , nv 



(18) 



are responsible for the effects of viscosity. An implementation of 
the diffusive fluxes in this conservation-law form is clearly the pre- 
ferred choice for FVM schemes, which are specifically designed for 
solving the integral form of these conservation laws. In fact, in this 
case they exactly conserve all the involved quantities to machine 
precision. We will therefore focus on this method in our study. The 
central aspect will be the numerical scheme used for estimating the 
velocity gradients at the cell interfaces, and hence the discretization 
of the diffusive fluxes. In the next section, we describe our approach 
for this in detail. 



3 A FINITE VOLUME SCHEME WITH VISCOUS 
FLUXES ON A VORONOI MESH 

3.1 Basic MUSCL-Hancock Finite Volume Scheme: 
Overview 

Finite volume methods enforce the integral form of the conserva- 
tion laws on discrete meshes. This approach is manifestly conser- 
vative, since fluxes of quantities that leave a cell simply enter the 



neighboring cell. The NS equations in finite- volume form are 

dQ* 



dt 



- 2_^ Aij Ff 
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dV 



(19) 



where, in general, the intercell fluxes contain both advective and 
diffusive contributions, 



■ adv, ij 



(20) 



The scheme used by AREPO is the finite volume MUSCL- 
Hancock approach, consisting of a MUSCL (Monotone Upstream- 
centered Schemes for Conservation Laws) linear reconstruction 
stage, and a Hancock two- stage time integration 



Qr +1 = Q?-At^ 






(21) 



hl/2 



where the numerical fluxes F^ represent appropriately time- 
averaged approximations to the true flux ¥ij across the interface 
shared by cells i and j. The time label n + 1/2 in Eq. |2T) in- 
dicates that an intermediate- stage (a half time-step evolution) has 
been performed to obtain the numerical estimate of F»j , meaning 
that the time-stepping in Eq. ( |2T) uses time-centered fluxes, giv- 
ing it second-order accuracy. The Hancock part of the scheme is 
a two-step approach (the familiar predictor-corrector algorithm) in 
which the correction half- step is obtained from the solution of the 
1-D Riemann problem across each face of the control volume. The 
general finite volume MUSCL-Hancock scheme has hence the fol- 
lowing three steps (Toro 2009): 

(I) Gradient Estimation, Linear Data Reconstruction and 
Boundary Value Extrapolation Once a local gradient estimate 
for the conserved quantities U* = (p, pv, pe)i of cell i is avail- 
able, linear data reconstruction takes the form 



jjR 



L U i +VU?(fi i -Si) 



(22) 



where we denote by U^ the estimated vector of conserved vari- 
ables at the centroid of the ij -interface, obtained by linearly ex- 
trapolating the cell-centered values U* of the i-th cell (on the "left" 
side) from s;, the cell's center position, to %. Similarly, U^ cor- 
responds to the estimates of the face-centroid values obtained by 
linear extrapolation of the cell-centered values of the j-th cell (the 
"right" side), whose center position is Sj. In both cases, Uj — fji is 
the position vector of the face centroid between the cells. The Jaco- 
bian VU™ is explicitly labeled with superscript n to point out that 
it corresponds to the estimate of spatial derivatives at the beginning 
of the time- step. 

(II) Evolution of Boundary Extrapolated Values This is, strictly 
speaking, the "predictor" half time-step. The conserved variables 
are evolved for At/ 2 with flux estimates obtained from the values 
at the beginning of the time- step: 

At 1 



U • = U 

13 lJ 2 V % 



u£ = uj 



At 1 



3 



F n 



(23) 
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(III) Solution of 1-D Riemann Problems and Computation of 
Godunov Fluxes This corresponds to the "corrector" half time- 
step in the two-stage Hancock approach. Once the values to the 
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right and left of the interface at time At/2 are known, the discon- 
tinuity is treated as a one-dimensional Riemann problem. An exact 
or approximate Riemann solver is used to return values of p, pw 
and pe at the interface, at a time corresponding to n + 1/2. From 
these values, the advective fluxes can be directly computed (Eq.[3]). 
These are time-centered fluxes F™ / used to update the system 
from the beginning of the time-step to its end, 



u- +i = ur,-Ati^^F- +i 



/2 



(24) 



Figures [T] and [2] illustrate the mesh geometry and the basic 
steps of this inviscid numerical scheme implemented in AREPO. 
One additional point we have not explicitly discussed here for sim- 
plicity is the treatment of the mesh motion, as indicated in Fig. [2] 
This is incorporated into the scheme by evaluating all fluxes in 
the rest frame of the corresponding face, as described by Springel 
(2010). This requires appropriate boosts of the fluid states and the 
fluxes from the lab frame to the rest frame of each face, and back. 
For a Voronoi mesh, the face velocities are fully specified by the 
velocities of all the mesh generating points. The latter can be cho- 
sen freely in principle, but if they are set equal to the fluid velocities 
of the corresponding cells, a Lagrangian behavior and a manifestly 
Galilean-invariant discretization scheme is obtained in which the 
truncation error does not depend on the bulk velocity of the sys- 
tem. 



3.2 A MUSCL-Hancock Finite-Volume Scheme with Viscous 
Terms 

A cell-centered, finite-volume solution of the NS equation can be 
written as 



Q n +1 = Q „ _ M £ ^p«+i/2 _ At £ Aij F2£(?, (25) 



where we have retained the distinction between advective and vis- 
cous fluxes. As in the case of the Euler equations, the numerical 
method essentially consists of the problem of finding accurate time- 
centered numerical fluxes across each of the interfaces of a cell. 
How to do this in detail for the diffusive part of the fluxes has been 
the focus of numerous efficiency and stability analyses (see Puigt 
|et al.|2010[ fo r a detailed description). 

Eq. (|25|) uses time-centered fluxes, obtained here with the two- 
step Hancock technique, as described above. Thus, for estimating 



^n + l/2 



^n + 1/2 



both F a J v (. and F di ^ (. a half time- step predictor stage is re- 
quired. In the MUSCL-Hancock approach for inviscid flow, this 
step is carried out by linear reconstruction from each cell center 
to the interface, followed by solving a one-dimensional Riemann 
problem at the interface where the extrapolations meet. The tra- 
ditional formulation of the Riemann problem and its solution are 
exclusive to hyperbolic differential equations and thus do not pro- 
vide exact solutions for the NS equations. Since a general solution 
for the viscous Riemann problem does not exist, we will treat the 
viscous fluxes in Eq. ( |25] > as a correction to the solution of an oth- 
erwise inviscid flow. 

Our NS version of the MUSCL-Hancock scheme consists of 
the following three different stages (in addition to those described 
in Section. [TT}: 

(A) Correct the MUSCL linear extrapolation of primitive variables 
by applying a viscous kick. 

(B) Extrapolate the cell-centered gradients linearly and evolve 
them for half a time- step. 



(C) Average the extrapolated velocity gradients at the interface 
and use them to estimate viscous fluxes. 

To extrapolate the gradients from their cell-centered values to 
the interfaces, information about the higher-order derivatives of the 
primitive variables is needed. If gradients are assumed to vary lin- 
early in space, an estimator for the Hessian matrix for each of the 
five primitive variables is sufficient. Evidently, enough information 
is contained in the cell-centered quantities to estimate both the local 
gradient V0 and the Hessian H^ corresponding to a given scalar 
quantity 0. However, estimating both of these simultaneously is 
significantly more difficult than estimating them one after the other. 
Therefore, we will effectively treat (j) and V0 as two independent 
fields that vary linearly in space, and this variation needs to be esti- 
mated from the mesh data through a suitably discretized differential 
operator. 

As a simpler alternative to the gradient reconstruction ap- 
proach, we briefly describe how one can use the gradients already 
available from the linear reconstruction step. In this approxima- 
tion, a given quantity varies only linearly within the control volume, 
such that consistently evaluated gradients are piece-wise constant. 
This means that each interface represents a discontinuity in the gra- 
dient field V0. Naively, one may think that the arithmetic average 
of both gradients that meet at a face is a good estimate for the gradi- 
ent at the interface itself. However, on second thought, one realizes 
that both cells do no necessarily have the same weight if cells of 
different volume meet. Furthermore, the unweighted average of the 
two cell-centered values really represents the value at the midpoint 
of the two mesh-generating points, which, for a Voronoi mesh, can 
be substantially offset from the mid-point of the face. We there- 
fore adopt the approach of Loh (2007), which consists in choosing 
one of the two gradients that meet at the interface, based on prior 
knowledge of the direction of the flow across the interface. Thus 
the three- stage scheme introduced above could be alternatively re- 
placed by the simpler method: 

(A- C) At the cell interface where two different gradients meet, 
choose the upwind gradient. 

In either method, once we have an estimate of both viscous and ad- 
vective fluxes, the time-step evolution of the conserved quantities 
Qi is carried out as in Eq. ( [25] ). However, the approach (A-C) is 
preferable to the Loh (2007) scheme because it uses time-centered 
estimates for both F™dv (. and F diff (. , hence preserving the or- 
der of accuracy of the original inviscid scheme. We therefore now 
provide a more detailed description of the individual steps in this 
three-stage approach. 



(A) Viscosity Kicks 

Although Eq. ( [25] ) is written in an unsplit form, the predictor step 
is indeed operator split, evolving the advective and diffusive terms 
separately (e.g. Coirier & Powell 1996). While our method for esti- 
mating the advective fluxes remains the MUSCL-Hancock scheme, 
the technique for estimating the diffusive fluxes is essentially con- 
tained in the estimation of the velocity gradients at each interface 
(see Coirier 1994, Puigt et al. 2010 for a series of tests on differ- 
ent interface gradient estimates). Looking for better accuracy, we 
have chosen to couple these two otherwise independent procedures 
by correcting/biasing the linear extrapolation of the velocity field 
(stage (/) in Section [XT) with a viscous source term. 

The benefit of carrying out a linear extrapolation to cell inter- 
faces in primitive variables is the simplicity of the Galilean trans- 
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Figure 1. Schematic representation of the mesh geometry and the MUSCL-Hancock integration scheme implemented in AREPO: a) TheVoronoi mesh is 
uniquely determined by the location of the mesh-generating points, b) A gradient estimate for all primitive variables is obtained from the immediate neighbors 
of a given cell, c) The gradient-estimation process is repeated for each cell in the domain and thus a piece-wise linear reconstruction is obtained for each 
primitive variable, d) The primitive variables are extrapolated toward each interface and evolved for half a time-step, e) For each face, a pair of extrapolated 
quantities for two neighboring cells i and j forms a local Riemann problem, f) The Riemann problem is solved for each face of a cell, yielding time-centered 
Godunov fluxes for the entire boundary of the control volume Vi of cell i. These fluxes are used for updating the conserved quantities of the cell through 
Eq.gl). 



formation needed to boost the quantities to the frame of a moving 
interface. Since the Galilean boost does not affect the mass and 
pressure of a given cell, only the local velocity field is transformed. 
In addition, adding force source terms to the equations of motion in 
primitive variable formulation is simpler, since these only couple 
to the momentum equations. Thus, a "viscous kick" can be applied 
to the velocity field in the half time-step evolution stage: 



Av visc = 



At 



^V 2 v + ^±i^V(V-v) 
LP P 



(26) 



In this way, the subsequent linear extrapolation of primitive vari- 
ables will already include viscosity effects to first order in time. 

While working with numerical fluxes across interfaces re- 
quires velocity gradients, the use of cell-centered source terms in 
Eq. \26\ calls for second order derivatives of the velocity field. 
Thus, in addition to the cell-centered velocity gradients \7v x , Vv y 
and Vv z , the cell-centered Hessian matrices H Vx , H Vx and H Vx 
are now needed. As we will see below, these matrices will be of 
use in more than one occasion, justifying the computational cost 
incurred to calculate them. 



(B) Linear Extrapolation of Gradients 

The linear reconstruction implemented in our MUSCL-Hancock 
approach essentially assumes that the gradient of a scalar quan- 
tity 4> does not vary significantly across the spatial scale of a cell. 
For smooth flows, the gradients of two neighboring cells V0 and 

V0 will not differ significantly. Furthermore, in the presence of 

strong discontinuities, gradients on each side will be slope-limited, 
and therefore will not differ by much either. Hence, a first guess for 
the gradient at the interface between two cells is just the average of 
the cell-centered estimates at each side of the face 



V0 = 



(V0)i + <v<^ 



(27) 



However, as we pointed out earlier, the gradient average above is 
actually representative of the midpoint between the two cell cen- 
ters Yi and Tj, which in general does not lie close to the center of 
the face in a Voronoi mesh, and may in fact lie within a third cell. 
Unless gradients are assumed to vary within a cell, it will not be 
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Figure 2. Detailed description of the flux calculation with a Riemann solver in step e) of Fig.^ For the case of a moving mesh, the standard MUSCL-Hancock 
method needs to be augmented with Galilean-boosts, as described by Springel (2010): (1) The extrapolation towards each interface is followed by a Galilean 
boost of the velocities to the rest frame of the face, and by a rotation of the coordinate axes. Each face is then treated as a one-dimensional discontinuity. Thus, 
the axes are oriented in the rotated frame such that the a/ -axis coincides with the normal to the face (left panel). (2) The primitive variables in the moving 
frame are evolved for half a time-step, including source terms if present (e.g. gravity or viscosity). (3) A one dimensional Riemann problem is solved at the 
interface. (4) The velocities are translated back to the lab frame and the advective fluxes are computed. 



possible to assign the estimate to the center of the interface with 
any confidence. 

Let us assume that the scalar field <fi(r) is infinitely differen- 
tiable and, consequently, so is its first derivative. Thus, we can Tay- 
lor expand both quantities to arbitrary order around a mesh gener- 
ating point 1*0 : 



M + ' 



(r-ro) 



+ i(r-r ) T H*| (r - r ) + 0(d 3 ) 

A lr 

V0(r)=vJ +H | (r-r ) 

lr lr 

+ J(r-r ) T D^| (r-r ) + O(d 3 ) 

A lr 



(28) 



(29) 



where H^ is the Hessian matrix of the scalar quantity <j> and D^ 
is a 3 x 3 x 3 tensor containing the third-order derivatives of (j) 
(i.e. D a bc — d 3 4>/dx a dxbdx c ). Truncating both Taylor expan- 
sions to first-order in d = r — ro, we see that we can obtain linear 
reconstructions for both the physical quantities and their gradients 
provided that we have numerical estimates for both the gradients 
and the Hessians at each mesh generating point. 

We emphasize that a Taylor expansion is not equivalent to a 
polynomial data reconstruction. Indeed, it is desirable that recon- 
struction schemes are manifestly conservative, in the sense that the 
average of the reconstruction over the cell should be identical to 
the value of 4> at the geometric center of the cell. This property of 
reconstruction schemes is sometimes referred to as if -exactness, 
meaning that if a polynomial reconstruction is cell-averaged over 
the mesh, the reconstruction procedure recovers the same polyno- 
mial. This condition is trivially satisfied for a linear reconstruction 
of the form 0(r) = & + (V</>)i(r — so). However, higher-order 
reconstruction schemes require the use of zero-mean polynomials, 



which, beyond first-order, differ from the Taylor series (e.g. Colella 
|& Woodward|1984||CoTrier & Powell|1996| ). 

The linear reconstruction of the scalar field (j) and of the vector 
field V0, treated as if they were independent quantities, effectively 
constitutes a hybrid method between standard linear reconstruction 
and fully K-exact second-order reconstruction, as illustrated in Fig- 
ure|3] In this approximation, second derivatives are considered neg- 
ligible for the spatial reconstruction of the primitive quantities, but 
they are still included for a more accurate estimate of the gradients 
near the cell interfaces. We also note, that in this way our numerical 
scheme reduces to that originally in AREPO (which is second-order- 
accurate) when the viscous fluxes are disabled. 

Once an estimate for the Hessian matrix H^ is available 

(Section [33] >, a linear extrapolation of the gradients from the cell 
centers to the interfaces can be obtained from 



(V0>i + <H*>(fy-r<), 



(30) 



which is a better approximation than Eq. ( [27] ). However, the time 
evolution of the gradients during a single step could be equally im- 
portant as their spatial variation over the length scale of a cell, 
hence we also need to evolve them for half a time- step to obtain 
a time integration scheme that is consistent with the second-order 
accurate two- stage MUSCL-Hancock approach. In the latter, to ex- 
trapolate and evolve a scalar quantity (j) we consider 

tj>\ =<t> l +V<t>\ (fy-sO-^/^ 01) 

\ij I r A \Ot I i 

where the time derivative of the quantity <j) in the control volume of 
the z-th cell can be obtained from the primitive variable formulation 
of the Euler equations in tensor notation: 

d t W a + A apb (W)d b Wp = 0. (32) 

Here sums over repeated indices are implied. Latin indices 
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Figure 3. Schematic representation of the double linear reconstruction proposed in this work compared to standard linear reconstruction and parabolic recon- 
struction. 



a, 6, c, cL. take the values 1,2,3 or x,y,z, while Greek in- 
dices a, /5, 7, ... take the values 1, 2, 3, 4, 5 and are used to num- 
ber the components of the primitive quantity vector (W a = 
p, v x , v y , v z , P for a — 1,2,3,4,5, respectively). As with our pre- 
vious notation, the indices z, j and k are reserved for labeling the 
mesh generating points and their associated cells. 

Eq. ( [32| is an advection equation for the primitive variables. 
Analogously, to "advect" the gradients of the primitive variables 
from the cell center to the interface, we can ignore the viscous terms 
and derive an equation of motion for the spatial derivatives by dif- 
ferentiating Eq. ( [32] ): 

OaOtWa + (d a A a p b ) d b Wp + AafrdadbWp = 0, (33) 

where we can identify the Jacobian matrix of the primitive variables 
as J oca = d a Wa — W a ,a, and the Hessian tensor (5 x 3 x 3) of 
the primitive variables as Hpba = dbdaWp = W/3,b, a - Therefore, 
the time derivative of each component of the primitive variable Ja- 



cobian matrix is 



dtJaa — B a (3t, a J(3b — A a (3bH{3 



(34) 



where we introduced the rank-4 tensor B a pba = d a A a pb — 
A a /3b,a- Since A a pb is a function of the primitive variables W a , 
the tensor B af 3ba can also be written as (see Appendix) 



B a pba — ^ d a Wry, 



(35) 



and therefore its numerical estimate is given by the product of the 
exact derivatives dA a pb/dW 1 (evaluated with values of the prim- 
itive variables at the center of the cell) and the (already available) 
numerical estimates for the gradients d a W 1 — J ia - The second 
term on the right hand side of Eq. (341 is the product of the known 
coefficients A a pb (evaluated at the center of the cell) and the nu- 
merical estimates of the Hessian tensor Hpba- 

Finally, with a numerical estimate of Hpba at hand (see Sec- 
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Figure 4. Sketch illustrating the individual steps involved in the extrapola- 
tion and half time-step evolution of the gradients, analogous to the advective 
flux calculation shown in Fig. [2] The different steps are: (1) spatial extrapo- 
lation of the gradients, followed by (2) a time advance by At/2 according 
to Eq. pi) , and (3) an approximate evaluation right at the interface. In step 
(4), the viscous fluxes are determined by evaluating Eq. fl8) with the values 
of the primitive variables and the velocity gradients at the interface. 



tion |3.3| ), the extrapolated and half time- step evolved gradients of 
the velocity are (in analogy to Eq.[3T): 

At / d Vv x ' 
~2~ 



V«J =(Vv x } i + (H Vx ) i {f ij 



+ 



dt 



(36) 



with analogous expressions for Vv y \ij and Vv z \ij. In Eq. ( [36] ), 
the term (dVv x /dt) i is obtained from Eq. (34i with a = 2 and 
a =1,2, 3. 

In Fig. [4] we show a sketch of the different steps involved in 
obtaining time-centered diffusive fluxes. We point out that taking 
the Hessian matrices of the velocity field to be identically zero is 
not equivalent to the alternative scheme {A'). The third term to 
the right hand side of Eq. ( |36| > is still different from zero even if 
Hpba = (Eq.[34j since, in general, B a /3b a J/3b / 0. By advecting 
the gradients according to Eq. (34) we gain additional accuracy at 
no additional computational expense because the terms B af 3baJ/3b 
are known exactly (see Appendix), given the values of the primitive 
variables and their respective gradients at the center of each cell. 

(C) Viscous Flux Calculation 

An accurate estimate of the viscous fluxes between two cells re- 
quires an accurate estimate of the velocity gradients at the inter- 
face. The gradient extrapolation method described above produces 
in general two different values of the velocity gradient that meet 
at the interface. This defines a general Riemann problem for the 
differential equation in Eq. d34j> which is no longer a homogeneous 
hyperbolic differential equation. Therefore, attempting to solve this 
new Riemann problem for the spatial derivatives of the scalar quan- 
tities introduces a significant additional difficulty. For simplicity, 



we will assume that the differences between two gradient extrapo- 
lations meeting at an interface are small enough such that a simple 
arithmetic mean can be used. This assumption, of course, is valid 
only when the field of second derivatives is sufficiently smooth (see 
Section[33]>. 

The time and area averaged flux across the face i-j that moves 
with speed w is defined as 



\l A. l3 J At J A .. L 



At 



■ 3idv,ij 



-F diff (W,<9W/<9r) dAij dt 



(37) 



Fdiff,2_ 



The numerical or Godunov estimate of these fluxes is chosen so 
that the analytic expressions for F a dv(U) and Fdifr (W, dW/dr) 
are evaluated with numerical estimates of U, W and dW/dr at 
the centroid of the interface. The advective Godunov fluxes are 



iv,ij — ■ 



F a dv(U Riem ) 



TT lab T 

U Riem w 



] ™i 



(38) 



where U Riem is the conserved variable vector at the centroid of the 
interface, as seen in the lab frame, obtained from the solution of a 
1-D Riemann problem across the i-j interface and along its normal. 
Multiplying by hij is equivalent to projecting the flux matrix F a d v 
(Eq.pl along the normal of each face. The Godunov fluxes F a dv, ij 
and Fdiff,ij are thus 5-component vectors. The diffusive Godunov 
flux vector is obtained from the diffusive flux 5x3 matrix 



v x Il x 





















llcccc 
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±±xy 
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Tlzx 
Kzy 

u zz 

V X U ZX + VyU Z y + V Z U Z 



(39) 

where U a b are the components of the viscous stress tensor II, 
which depend on the local value of the velocity and the velocity 
gradients. These components are: 



J-J-xx 


= 


-r]d x v x - -n (dyVy + d z v z ) + ("V • v 


i^yy 


= 


4 2 

-ndyVy - -n(d z v z + d x v x ) + CV • v 


u zz 


= 


4 2 

-r)d z v z - -n(d x v x + dyVy) + (V • v 

6 6 (40) 


J-J-xt/ 


= 


n yx = rj (d y v x + d x v y ) 


n„* 


= 


Hzy = r\(d z v y + d y v z ) 


u zx 


= 


Rxz = r\ (d x v z + d z v x ) 


Just like with the advective fluxes, the flux tensor (Eq. [391 must 
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be projected onto the normal hij of each ij -interface to obtain the 
5 -component vector 



because 



= F d 



Wkt m ,(^W/ar)i a D b pi 



<)nii, 



(41) 



where W^m is the primitive variable vector at the centroid 
of the interface, as seen in the lab frame (whose associated 
conserved variables are UrL in Eq. [38b. The spatial deriva- 



\lab 



tives (dW/<9r)^ D x correspond to our extrapolate-and-average 
scheme for linearly varying gradients. As with W^, we are in- 
terested in estimates of dW/dr at the centroid of the face. For 
both these quantities, only the velocity and its spatial derivatives 
are relevant when viscous fluxes are calculated. 



3.3 Hessian Estimation 

In analogy to the gradient calculation for Voronoi meshes dis- 
cussed by Springel (2010), here we discuss the estimates of the 
cell-centered Hessian matrices for each of the primitive variables 
W a . To this end, let us consider a vector field u that varies approx- 
imately linearly with distance as u « Ui + h (r — n) near r*. Up 
to linear order, the first derivative of u is simply h. The volume- 
average of the spatial derivatives of u in the vicinity of r$ is 



V t 



du\ 



du 



dV 



L 

Jd\„ 

3& Ai i 



fVi ® r 
udA 



(42) 



+ h(r - n 



-cL4, 



where we have assumed that the linear approximation is valid up to 
all the neighboring mesh generating points Yj . It is straightforward 
to verify that the average matrix (du/dr)i can be written as 



du 
dr 






Ui + u ? 



vE^' ( hc 



V, 



jy^i 



(43) 



Writing the vector product ( A u) v in tensor form (where A 
is a n x n square matrix and u and v are vectors of dimen- 
sion n), it is easy to prove the identity A ac u c Vb = A ac v c Ub + 
£bfc£fdeUdV e A ac . Equivalently, going back to vector notation, we 
have ( A u) ® v = ( A v) u + (u x v) x A, where, for simplic- 
ity, we used vector notation to denote a "cross product" between a 
vector and a matrix. 

Therefore, the second term on the right hand side of Eq. ( [43] ) 
can be written as 



y^Aij (hcij <g> ^ ) =V^ (hr 



+£M 



J^i 



(44) 



Here, the second term on the right hand side vanishes identically, 






J\i- 



x h 



{/ (• 



Vx 



ri 


+ r A „ 




2 )" 


r - 


Ti + Yj 



x dA xh 



V) 



dy xh 



= 



(45) 

On the other hand, the first term on the right hand side of Eq. \AA) 
can be rewritten by means of the replacement h r^ = — h (yj — 
ri) = Ui — Uj. Finally, identifying the vector u; with the gradient 
(V0)i of a scalar quantity 0, and the matrix (du/dr)i with the 
cell-averaged Hessian matrix (H^)i, Eq. J441) takes the form 



w 



3^i 



+ Wh 



/ Vij ) 



(46) 



The most noteworthy characteristic of this expression is that it is 
purely algebraic and explicit in nature. That is, the Hessian matrix 
of (j) is simply a linear combination of the neighboring gradients 
in which the coefficients are predetermined quantities that depend 
only on the local mesh geometry. Each one of those neighboring 
gradients is, at the same time, a linear combination of its imme- 
diate neighbors' scalar quantities (see Eq. 21 of Springel 2010). 
Therefore, the Hessian estimate of Eq. ( |46| is a weighted linear 
combination of scalars from its immediate neighbors and from its 
neighbors' neighbors and, as such, it implicitly employs a larger 
stencil than the one used for the gradients. 



3.4 Slope-Limiting the Hessians 

It is well known that higher-order reconstruction schemes are prone 
to produce spurious oscillations in the vicinity of steep gradients, 
unless this is prevented by appropriate slope limiter methods (Toro 
2009 ). These non-linear corrections in the data reconstruction pro- 
cedure prevent overshoots and allow for true discontinuities in the 
solutions. So far, we have discussed how to estimate second deriva- 
tives from first derivatives, which in turn are already reconstruction 
estimates obtained from the scalar physical quantities. Potential ir- 
regularities in the second derivative fields can lead to spurious os- 
cillations and unphysical values of the viscous stress tensor at the 
cell boundaries. To alleviate this problem, we enforce local mono- 
tonicity of each component of the gradients, which is equivalent to 
smoothing out the Hessian estimates. In practice, this is achieved 
by replacing the Hessian matrix by a 'slope limited' version 



W ; =A,H 



(47) 



with 



OL X i 











a y 











al 



and where o% — min (1,-0^ 

(48) 
are the slope limiters < af < 1 for each direction x, y and 
z. This MINMOD-type slope-limiting method is readily applicable 
for irregular meshes. The quantities ipij are defined as 

f ({d a <t>)? !iX -{d a 4>)i)/A(d a (i > )^ A(d a 4>) 



^ 



((dafii 



(da^A /A (0„0. 



A(da4>) x 

A(^) 2 



>0 
<0 
= 
(49) 
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where the A (d a 0) i7 - are the components of the vector A (V</>) • ■ = 
(H (f) )i(fij — Si), i.e. the estimated change in the gradient V0 be- 
tween the centroid f^ of the cell and the center of cell i. The quan- 
tities (<9 a 0)™ ax and (<9 a (/>)™ m are the maximum and minimum of 
the a-th component of the cell-centered gradient estimates among 
all neighboring cells of cell i, including i itself. 

To our knowledge, the slope-limiting technique has not been 
applied to the second derivatives before. However, its purpose is 
equivalent to the "flattening" procedure near shocks carried out by 
|Colella & Woodward] < [1984| ) for parabolic reconstruction. In our 
approach, the suppression of oscillations near shocks is exclusively 
handled by the limitation of the gradients, since the reconstruction 
of hydrodynamic quantities is only of linear order. Thus, the Hes- 
sian limitation procedure serves the sole purpose of guaranteeing 
a smooth variation of the gradients and avoiding spuriously large 
viscous fluxes. Future improvements of the present method could 
however also employ these second derivatives for higher-order re- 
constructions of the scalars. 



3.5 Time integration and time-step criterion 

Because of the more complex mathematical properties of the NS 
equations compared with the Euler equations, obtaining a rigorous 
analytic expression analogous to the CFL stability criterion for the 
allowed time step size is not possible. However, MacCormack & 
|Baldwin| fl975) obtained an approximate semi-empirical stability 
criterion when advective, viscous and heat diffusion terms are con- 
sidered. When there is no heat flux, the time- step criterion can be 
written as (e.g. |Kundu & Cohen|2008] > 



At < 



crAt 



CFL 



1 + 2/{Re}, ' 



(50) 



where A£cfl is the standard CFL-criterion time-step except for 
the Courant-Friedrichs-Levy coefficient, which is absorbed into a 
"safety factor" a (usually w 0.9). In Eq. ( [50| , the cell Reynolds 
number {Re}* is 



{Re}* 



pWi\Ri 



(51) 



where v^ is the velocity of the gas relative to the motion of the 
grid and Ri is the effective radius of the cell, calculated as Ri — 
(3Vi/47r) 1/s from the volume of a cell (or as Ri = (Ai/n) 1/2 
from the area in 2D). Similar approaches to derive an appropriate 
NS time- step have also been described by Mavriplis & Jameson 
( fT990] > and |Coirier & Powell|fl996l . 

The numerical integration scheme we employ is time unsplit, 
that is, advective and diffusive fluxes are applied simultaneously 
during each hydrodynamic time-step and not sequentially (Eq.[25]>. 
The prediction stage, on the other hand, is operator- split, since the 
advective and diffusive terms are computed almost independently 
of each other. This is in part due to the nature of the standard one- 
dimensional Riemann problem, whose solutions - strictly speaking 
- are only valid for the hyperbolic Euler problem, but are not solu- 
tions to the full NS equations with their additional parabolic terms. 
The validity of this approach ultimately relies on the assumption 
that the viscous terms in the NS equations are typically small per- 
turbations to the Euler equations. 



4 NUMERICAL TEST RESULTS 

To test the performance of AREPO when our new treatment of vis- 
cous fluxes is included, we have carried out a number of test sim- 
ulations for physical situations with known analytic or quantitative 
solutions. Usually, the problems with known exact solutions are 
either of self- similar type or have symmetries that make the non- 
linear term proportional to (v • V)v vanish identically. Owing to 
these limitations, numerical simulations of situations with exper- 
imentally well-established behavior, such as flow past a circular 
cylinder, have become common-place in testing the performance of 
NS codes. We will therefore also carry out such qualitative bench- 
marks, besides looking at a few simple problems with analytic so- 
lutions. 



4.1 Diffusion of a Vortex Sheet 

A simple problem of laminar flow in the presence of viscosity is 
given by the vortex sheet diffusion test. In this problem, the initial 
velocity field at t — is given by v = (it, 0, 0) with u — 1 for 
y > and u — — 1 for y < 0. Because of the symmetry of the 
problem, the NS equations reduce to a ID diffusion equation 



du 



<9V 



with solution (e.g. |Kundu & Cohen|2 008 ) 



erf 



2y/Hi. 



du e~ y l* vt 



dy 



(52) 



(53) 



In Figure [5] we show the time evolution we obtain for a 
two-dimensional simulation domain with initially uniform pressure 
and density (p = P = 1), and with a velocity field given by 
v = (sgn(?/), 0, 0). The mesh generating points were distributed 
regularly at the initial time to produce a Cartesian mesh. As the 
system evolves, the velocity and the vorticity fields as a function 
of time and vertical coordinate y follow the exact solution remark- 
ably well. It is worth pointing out that the initial singularity in the 
vorticity field is unresolved numerically (and thus appears as being 
uniformly zero throughout the domain), since the system is started 
with an exact sharp discontinuity. Static, perfectly aligned meshes 
with slope limitation techniques will typically maintain this unre- 
solved vorticity and thus no diffusion will proceed unless some nu- 
merical perturbations are seeded that break the mesh alignment of 
the initial state (a common way to overcome this difficulty is to 
start the system according to Eq. ( [53] ) at t > such that there is 
initial vorticity). However, the moving mesh of AREPO "sees" a 
non-zero velocity gradient as soon as the upper and lower halves 
of the domain become unaligned with respect to each other. This 
happens because, as soon as a cell shifts its position, the number of 
its neighbors that have a drastically different velocity increases and 
so does the "statistical weight" of the discontinuity. At this point, 
the slope-limiting technique, which had ignored the discontinuity 
in the perfectly aligned mesh, now identifies the local variation as 
"real" and the vorticity field is "detected". 

Fig.[6]shows the corresponding two-dimensional velocity field 
of the diffusing vortex sheet test at four different times, together 
with the geometry of the underlying Voronoi mesh. The mesh ge- 
ometry nicely shows how the cells transform from a Cartesian con- 
figuration to an unstructured mesh, while the velocity field evolves 
from a piece-wise constant state with a central discontinuity to a 
smoothly varying shear flow due to the effects of viscosity. 
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Figure 5. Diffusion of a vortex sheet. The two panels show the velocity u along the x-axis {left panel), and the vorticity {right panel), at times t = 0, 0.1, 0.2, 
0.4, 0.8, 1.6 and 3.2 (from black to red), for a dynamic viscosity coefficient /i = up = 0.005. The solid lines are given by the analytic solution described 
by Eqs. J53j, while the solid circles are all 2500 cell-centered velocity and vorticity values of the initially Cartesian 50 x 50 mesh. Note that the simulation 
is started with a sharp discontinuity in velocity and thus the 5-function vorticity field is initially unresolved. If the mesh would remain exactly Cartesian, the 
diffusion of vorticity would actually be suppressed in this case. Nevertheless, the small asymmetries introduced by the moving mesh trigger the diffusion 
regardless of the initially unresolved setup, and the time-dependent numerical result closely follows the expected exact solution. 



(54) 



4.2 Diffusion of a Gaussian Vortex 

The two-dimensional circular velocity distribution corresponding 
to an irrotational vortex of circulation F is 

where the vorticity uj = |V x v| = (l/R)d(Rvo)/dR is zero 
everywhere except at the origin (u = S(R), i.e. a vortex line). In 
a viscous fluid, this velocity profile has to be sustained by a point 
source of vorticity at the origin (e.g. an infinitely thin rotating cylin- 
der) otherwise the vortex line will decay in a similar way as the 
vortex sheet in the previous example. If the velocity at the origin is 
set impulsively to zero, the subsequent evolution of the azimuthal 
velocity is given by 



v e (R,t) = 
while the vorticity uo — Y 



-R z /4 



2ttR I 
V x (v 6)\ 



vt\ 



z evolves as 



-R z /4ut 



(55) 



(56) 



6 



(57) 



47VVt 

and the Laplacian of the velocity field is 

|Vv| = 2^(2^F e 

Because of its geometry, this problem is significantly more chal- 
lenging than the vortex sheet test considered above and cannot be 
impulsively started at precisely t = 0. Besides the initial singularity 
in the vorticity field, the velocity field is divergent as we approach 
the origin. In addition, it is not possible to capture the azimuthal 
velocity field when the distance from the origin is comparable to 
the grid resolution. At the same time, the azimuthal velocity field 
is challenging for the boundary conditions, because the problem is 



self-similar in nature and therefore natural boundaries do not exist. 
These problems did not exist for the vortex- sheet problem, which is 
of one-dimensional nature. Nevertheless, evolving the system from 
an initial time t > minimizes most of these complications. In 
addition, we extend the computational domain far beyond the re- 
gion of interest, such that boundaries become essentially irrelevant 
during the timespan of the numerical solution. 

We setup a Cartesian mesh (100 x 100) with an imposed initial 
velocity profile of 

I. — exp 



ve,o 



r 

2^R 



with v ■ 






(58) 



r z y 

4zyt 7 _ 

corresponding to a Gaussian vortex that we center in the middle 
of the domain, which extends over the range [0, 40] x [0, 40], and 
thus accommodates a radial range from R = to R = 20. The 
adopted physical parameters are to = 10, fi = 0.08, r = 1.0, and 
the initial density field is constant with p = 1. The pressure field, 
however, is not uniform because the fluid is not started from rest. 
We obtain the correct pressure profile from the radial component of 
the equation of motion: 

_v\ _ _ldP 
R ~ pdiT 
and thus the initial pressure profile is 

1 (^ ( R 2 \ 



-R 2 /(2ut ) \R 2 /(4ut ) 



1 



+ - 



Ei 



Ei 



4vto V~~ V 2zy ^o ) ^~ V 4 ^o 
where Po is an integration constant. The precise value of Pq is 
irrelevant for the similarity solution presented here, because it is 
obtained for incompressible flow. In our numerical experiments 
(which are compressible), we set Po such that P = 1 at R = 0. 
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Figure 6. Time evolution of the mesh geometry and the velocity field for a diffusing vortex sheet test. As the vorticity spreads from the center of the domain to 
the upper and lower boundaries, the mesh adapts to the continuous change in velocity until its original Cartesian structure disappears entirely. The color table 
(from blue to red) corresponds to the range between u = —1.0 and u = 1.0 in linear scale. 



Fig.|7]shows the time evolution of the velocity field, the vortic- 
ity field and the Laplacian field for a Gaussian vortex started on an 
initially Cartesian mesh. We find not only that the velocity evolves 
as expected based on the similarity solution, but the first and second 
derivatives also show excellent agreement with the analytic expec- 
tations. These results validate both the space- and time-accuracy of 
our viscous integration scheme, as well as the accuracy with which 
the second derivatives are estimated. 



4.3 Plane Poiseuille and Couette Flows 

Next, we consider impulsively- started plane Poiseuille and Cou- 
ette flows where a fluid between two parallel plates is initially 
at rest, and then, suddenly, either pressure gradients or plate mo- 



tions are applied. The time-dependent solution has the form v = 
(u(y,t), 0,0), where the horizontal velocity can be decomposed 
into steady and time-dependent parts, u(y, t) = m(y) + u(y, t). 
In the presence of a pressure gradient and an upper plate moving 
at constant speed U, the steady state solution is the well-known 
expression 



uo(y) = 



yU _ y_dp 

b 2/jL dx 



(p-y), 



(59) 



for which the special cases [7 = and dp/dx = are commonly 
known as plane Poiseuille flow and plane Couette flow, respec- 
tively. 

The time dependent component u(y, t) is a solution of 
Eq. (|52|, subject to the initial condition u(y,0) — —uo(y) and 
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Figure 7. Time evolution of a diffusing Gaussian vortex. For each time (as labeled), we show the azimuthal velocity profile v$ (R), the vorticity profile uj(R) 
and the Laplacian profile V 2 vq, as computed by AREPO (blue points; only a random 10% of the total shown) and compare it to the corresponding analytic 
expressions (solid red lines). 



the boundary conditions u — at y — and y — b. By separation 
of variables, the general solution is (e.g |Graebel|2007| ) 



u(y,t) = y^ A n e 



ut/b 2 



nTiy 



(60) 



where the coefficients A n are determined by the initial condition 
/>o(2/)sin^ch, 



A- n 



rb . 

Jo S1 



2 mvy 



dy 



2U(-l) n 2 dp 



bji dx \ nTT 



(61) 



[l-(-l) n ] • (62) 



The numerical setup for this problem is straightforward. We 
produce a Cartesian mesh in the range [0, 1] x [0, 1] with a resolu- 
tion of 50 x 50. The fluid is originally at rest and its density and 
pressure are given by p — P — 1. The equation of state is that of an 
ideal gas with adiabatic index 5/3. To represent the plates, the up- 
permost and lowermost rows of cells are replaced by "solid bound- 



aries" at which the no-slip condition is enforced, i.e. v x 







(see Fig. [10]). Moving solid boundaries are straightforward to im- 
plement with a Voronoi tessellation mesh. A solid surface can be 
constructed as a series of mesh-generating point pairs, one on each 
side of the surface, such that the common interface - equidistant to 
both points - defines the boundary locally (see Serrano & Espanol 
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Figure 8. Impulsively started plane Poiseuille and Couette flows as a function of time, a) Time evolution of the horizontal velocity profile versus vertical 
distance. Solid curves represent the analytic solutions of Eqs. |59J to (61) for U = and dp/dx = —0.05, at ten different times (time increasing from black 
to red). The data points correspond to all the cell-centered values of velocity along x for a 50 x 50 simulation started from rest, b) Time evolution of the 
horizontal velocity profile versus vertical distance. Solid curves represent the analytic solutions in Eqs. {59| to (61) for U = 0.1 and dp/dx = at eight 
different times (time increasing from black to red). The data points correspond to all the cell-centered values of velocity along x for a 50 x 50 numerical 
simulation started from rest. 



|2001| and|Springel 2010). The Voronoi cell on the side of the "solid" 
object can regarded as "a ghost cell within the domain". That is, 
this cell is part of the domain discretization process and is tessel- 
lated/updated as any other normal gas cell. However, when solving 
the Riemann problem at the local interface between a "solid" cell 
and a real gas cell, boundary conditions are imposed in the same 
way as boundary conditions on the outer box are imposed. For per- 
fectly reflecting boundaries, the normal component of the velocity 
is reflected in the "solid side" or "outside region" of the interface. 
For non-slip boundaries, the entire velocity vector is reflected, such 
that the velocity at the interface is zero (Figure [TO). 

We run two different test problems. For the first one, both 
plates remain at rest and an external gradient of dp/dx = —0.05 is 
imposed. For the second test, the bottom plate is at rest and the up- 
per plate moves at a constant speed of U = 0.1. In both test simu- 
lations, the dynamic viscosity coefficient has been set to /n = 0.05. 
In Figure [8] we show the time evolution of the horizontal veloc- 
ity profile both for the plane Poiseuille and Couette flows. In both 
cases, the numerical results match the analytic expectations very 
well. In Figure ^9] we also show maps of the velocity profile and the 
mesh geometry at different times for the Poiseuille case. The grid 
evolution shows how the Cartesian structure is progressively lost, 
but that the dynamic Voronoi mesh of AREPO successfully avoids 
any mesh-tangling effects. 



4.4 Time-Dependent Circular Couette Flow 

We now turn to a more challenging problem, which highlights the 
ability of our scheme to deal with geometrically complex boundary 
conditions. For purely azimuthal motion, the NS equations in the 



radial and tangential directions are 

_v\ _ _ldP 
~# ~ ~])dR 



dt 



= M 



dR 



]^_d_ 

RdR 



(Rve) 



(63a) 



(63b) 



The exact solution of steady flow (i.e. dve/dt = 0) between con- 
centric cylinders with boundary conditions vq = QiRi at R = R±, 
and v e = £l 2 R 2 at R = R 2 is given by (e.g |Kundu & Cohen|2008] > 

(VL 2 Rl - QiRf) R 2 - (H 2 - ft) RlRl 



v 9 ,o(R) = 



(64) 



R{Rl-R\) 

where Ri and ft (i = 1 , 2) are the radii and angular velocities of 
the respective cylinders. 

The impulsively- started version of this problem can be solved 
analytically by separation of variables (see Tranter 1968, Graebel 
2007). The full solution can thus be written as ve(R, t) = ve,o + 



vq (i?, £), where the time-dependent part has the form 

v d (R,t) = ^2{C 2 Ji(nR) + C 2 Y 1 (nR)}e- un2t , 

and J\ and Y\ are Bessel functions of the first and second kind, re- 
spectively. This time-dependent component is subject to the bound- 
ary conditions ve(R,t) = at Ri and R 2 , thus allowing us to 
eliminate C 2 : 



v e (R,t) = J2 



-B 1 (n s R)e- 



^Y^risR! 

where the n s are the roots of the equation B\ (nR) = with 
BiinR) = J 1 (nR)Y 1 (nR 1 ) - Y 1 (nR)J 1 (nR 1 ). 

Finally, the coefficients A s are determined by imposing the 
initial condition vq — —ve,o at t = 0. To solve for each coefficient 
independently, the steady state solution must be written in terms of 
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Figure 9. Time evolution of the mesh geometry and the velocity for flow between parallel plates. The horizontal velocity field u for plane Poiseuille flow is 
rendered at four different times. The evolution of the velocity field (see Fig 1 8 (a)) is accompanied by the evolution of the mesh from an initially Cartesian set 
up (top-left panel) to a fully unstructured grid by the time the flow has reached steady state (bottom-right panel). The (linear ) color scale ranges from blue 
(u = 0)tored(ii = 0.12). 



a series expansion of ve,o in the basis functions Bi(n s R). After 
some algebraic manipulations, we obtain 

Ji(n s R 2 ) 



V0 



,o(R) = 7TQ 2 R2J2 



' 1 J?(n a Ri)-J?(nsR2 



-B^risR) 



J\(n s R\) - Ji(n s R 2 



Qi Ri 



n 2 R 2 

and therefore the time-dependent component is given by 

~ fr> j.\ n V^ Ji(n s R 2 ) , . -un 2 t 

MRS) = ~^E Jl{nsRi) _ j ?{nsR2) Mn s R)e 

x [ VL 2 R 2 J 1 {n s Ri) - Ji{n s R 2 )£liRA . 



Collecting these results, the complete expression for the time- 
dependent angular velocity profile is 



Q(R,t) 



V0 TV 



£■ 



Ji(n s R 2 ) 



-B 1 (n s R)e- 



R RJ^JHnsRi)-JHnsR2) 

[ Q 2 R 2 J 1 (n s R 1 ) - J 1 {n s R 2 )VtiRi 



(65) 

We realize the moving boundary conditions in the present case 
through special Voronoi-cells with prescribed motion and boundary 
conditions, as described in Springel (2010). In the present case, we 
use two sets of mesh-generating points, each one consisting of a 
series of outside-inside pairs located on either side of the bound- 
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Figure 10. Schematic representation of (a) reflective and (b) non-slip boundaries within the computational domain. After the spatial and temporal extrapolation 
steps in the MUSCL-Hancock method (panel e) in Figure [TJ, the Riemann problem is solved as elsewhere in the domain but with the boundary-side cell 
mimicking the gas side with either one velocity component - the normal one - reversed (reflection) or all three (non-slip). 
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Figure 11. Geometry of the circular Couette flow. The left hand panel shows a schematic view of the two-dimensional problem. The right hand panel displays 
the actual initial mesh used in AREPO in a setup where we start the problem impulsively from rest. Each cylindrical boundary (at radii R\ and R2) is 
generated by two layers of cells, one side representing the solid cylinder and the other representing the fluid. These layers of cells are positioned along circles. 
The remainder of the fluid cells, originally at rest, are distributed like a Cartesian grid. The cells outside the outer cylinder are "auxiliary cells" and are only 
included to fill the computational domain, but do not exert any influence on the fluid inside the cylinders. The motion of the cylinders is prescribed to remain 
constant (with angular velocities fl\ and Q2X and thus represents a source of kinetic energy. The motion of the fluid in between the cylinders is induced by 
means of the no-slip boundary condition at the contact surface, and the momentum that is transported in the radial direction through the shear viscosity. 



ary and running parallel to it, so that two circular boundaries of 
radii Ri and R2 are defined which can be made to rotate at angular 
frequencies Qi and Q2, respectively. Note that the only significant 
technical difference between this problem and the preceding exam- 
ples is the way the boundary cells are prescribed to move; the rest 
of the numerical scheme remains unaltered. 

Figure[TT|illustrates the geometry of the circular Couette flow, 
and our realization of a suitable mesh in AREPO. Since the equa- 
tions of motion are always solved in the moving frame of the inter- 
faces, there is no practical difference between stationary and mov- 
ing boundaries when they are constructed as a part of the mesh. 



Figure [12] shows an enlargement of the mesh at the boundary cor- 
responding to the inner cylinder, which is represented by a set of 
Voronoi faces that follow a circular path. Each one of these Voronoi 
faces is defined by two mesh generating points located on either 
side of the face, one of them outside the cylinder on the fluid side, 
the other inside the cylinder on the side that does not contain fluid. 
The right-hand panel of Fig.[l2]shows the same region again, but at 
a slightly later time. This gives a sense of how the initial Cartesian 
mesh between the cylinders reacts to the fluid motion. Since the lat- 
ter is azimuthal, the mesh eventually develops an axial geometry, 
independent of the initially Cartesian setup. 
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Figure 12. A zoom showing the detailed mesh geometry around the inner 
boundary of the circular Couette flow, at two different times. The left hand 
panel shows a close-up view of the right panel of Fig. [TT] The Voronoi 
faces that make up the cylindrical boundary are created by close pairs of 
points, which either lie inside the solid cylinder (red) or on the fluid side 
(black). The gray cells define the contact region of the fluid domain with the 
cylinder; here the no-slip boundary conditions are imposed. An analogous 
geometry applies for the outer cylinder. The panel on the right hand side 
shows the same region of the computational domain at a slightly later time, 
when the mesh filling the fluid region has started to react to the motion of 
the cylinder. 



Our numerical experiment for this setup has the following pa- 
rameters. The initial mesh as described in Figs. [TT| and [T2] contains 
3,254 points, out of which 2,644 are regular fluid cells, 250 are 
boundary fluid cells, 250 are solid boundary cells and 110 are un- 
used auxiliary cells that are only put in to fill up the total mesh 
area to an enclosing rectangular shape, as presently required by 
AREPO. The radial distance between the cylinders is spanned by 
20 cells. The physical parameters of the Couette flow are Ri — 1, 
R2 = 2.5, Hi = 0.5, and Q2 = 0.1, with a dynamic viscosity 
coefficient set to \i — 0.005. In addition, since the flow is started 
from rest, the pressure and density are taken to be uniform with val- 
ues p — P — 1. Figure p3] shows the time evolution of the angular 
velocity profile as it asymptotically converges to the steady state 




Figure 13. Angular velocity profiles at different times for an impulsively 
started Taylor-Couette flow. For seven snapshots at times t = 0.5, 20, 40, 
60, 80, 100, and 120 we show the cell- centered values of Q, which are 
plotted as filled blue dots for all fluid cells in the calculation. No binning 
or averaging has been performed. The clustering of cell-center points as the 
system evolves is simply a consequence of the mesh adopting an axial sym- 
metry in an adaptive fashion. The dashed lines give the time-dependent an- 
alytic solution of Eq. {65) at the corresponding times. The numerical results 
are almost indistinguishable from the exact solution. The red curve depicts 
the steady state solution to which the time-dependent solution eventually 
converges. 



solution. The agreement of the numerical data points with the exact 
analytic solution (Eq.[65]> is exceptional at all times. 

Finally, we show in Fig. [14] the mesh geometry at the end of 
the calculation. Even though we have started the calculation with 
an initially Cartesian mesh, the memory of this geometry is lost 
during the calculation, and the mesh dynamically adapts to the az- 
imuthal flow structure present in this problem. The transition from 
a Cartesian grid towards a cylindrical-like mesh can also be seen in 
the output sequence of the simulation shown in Fig. [13] where the 
values of the radial position of the cells start to segregate into a set 
of radial "bins". The number of these radial clusters corresponds to 
the average number of cells along the radial direction. 

It is interesting to comment on the scatter of points - espe- 
cially at the beginning of the simulation - as seen in the angular 
velocity profile of Figure pj] This is a reflection of the challenging 
initial mesh geometry. Although high-order schemes - fifth or sixth 
order - are not sensitive to the compliance of the mesh geometry 
with the flow, second order schemes are. In this particular case, an 
axially symmetric mesh geometry would be more suitable due to 
the characteristics of the flow. However, the main point of this test 
is to show how the mesh responds to the evolution of the problem, 
achieving rough axisymmetry despite the unfavorable initial setup. 

As discussed by Springel (2010), our moving Voronoi mesh 
technique needs a "quality control" to keep cells sufficiently regular 
in order to avoid large errors in the spatial reconstruction. However, 
this modification of the mesh motion comes at a price: imagine a 
very strong compression along one direction (e.g. due to a very 
strong shock), then the mesh cells will acquire locally a high as- 
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Figure 14. Mesh geometry for the circular Couette flow towards the end of 
the numerical integration. Even though we have started the calculation with 
an initially Cartesian mesh, this structure is quickly lost in favor of an on 
average azimuthal mesh geometry. 



pect ratio, which our mesh-quality control motions will try to elim- 
inate, if needed by breaking the mesh symmetry (cell shapes are 
made "round" through small transverse motions). This is what hap- 
pens when we start the Couette flow impulsively on a non- suitable 
mesh. The introduction of asymmetries in the mesh can influence 
the flow, in particular in situations where fluid instabilities develop 
(see the Kelvin Helmholtz instability test in Springel 2010), where 
such asymmetric discretization errors can source growing pertur- 
bations. We note however that also on regular Cartesian meshes 
similar "grid-sourcing" errors exists. It appears unlikely that the 
poorer ability of the dynamic Voronoi mesh to maintain strict mesh 
symmetry is particularly detrimental for physical applications. 

4.5 Flow Past a Circular Cylinder 

We next consider the flow over a circular cylinder immersed in 
a wind tunnel. The geometric setup of the problem is shown in 
Fig. [15] The flow comes from the left at a fixed horizontal veloc- 
ity U. The upper and lower boundaries are also kept at constant 
velocity U. Far from the cylinder, at the right end of the compu- 
tational domain, we impose again an exit velocity U. The injec- 
tion and exit regions are forced to have the prescribed inflow and 
outflow mass fluxes at all times, something that we numerically 
impose through small "buffer" regions as labeled in Fig. [15] For 
static Cartesian grids, this buffer region does not need to extend 
more than one cell in the x-direction. However, moving grids re- 
quire not only the injection of momentum from the left, but also 
the injection of new mesh-generating points, since the wind tunnel 
will otherwise produce a depletion of cells at the left end as the 
mesh generating points drift to the right in the direction of the flow. 
We address this issue by letting the mesh automatically generate 
new cells through cell splitting, as already implemented in AREPO 
(see examples in Springel 2010). In doing this, some attention must 
however be paid to guarantee that the new cells reproduce the ex- 



ternally imposed inflow boundary conditions, which is most easily 
achieved with a sufficiently broad buffer region on the left end of 
the wind tunnel that covers the region where new cells are injected. 
Similarly, we employ the ability of AREPO to automatically remove 
mesh cells to prevent them from piling up on the right end of the 
wind tunnel. Altogether, we have created a wind tunnel that is filled 
with a mesh that blows with constant velocity from left to right, in 
a quasi- stationary state. 

The other geometric parameters of the test problem we sim- 
ulate here are the diameter d of the cylindrical obstacle, the width 
W of the tunnel and its length L. We have chosen W — 6.25 d 
and L = 5 W = 31.25d , and have scaled all length units such that 
W = 1.0. The flexibility of the Voronoi mesh allows us to easily 
embed a cylindrical obstacle within the initially Cartesian back- 
ground grid that fills the tunnel. Fig. [16] shows how we can tailor 
the mesh construction to reproduce the curved surface of the cylin- 
der, using techniques similar to those that we used for the circular 
Couette flow problem. 

The physical properties of the problem are primarily deter- 
mined by the external velocity of the flow, U, and the dynamic 
viscosity of the fluid p. In our numerical experiments we set the 
external flow velocity to U — 0.5, and combine this with constant 
initial pressure and density (p = P = 1). We take the fluid to 
be described by an ideal gas equation of state with adiabatic index 
7 = 5/3. The characteristic Reynolds number of the problem can 
then be defined by 



Re 



Ud 

v 



Udp 



(66) 



where p might however vary in time and space since the flow is 
fully compressible. 

We have performed several numerical experiments of this 
problem using the viscous module added to AREPO. In each of 
these simulations, the Reynolds number is the only relevant quan- 
tity being changed. This is accomplished by changing p exclu- 
sively, while keeping the other parameters fixed. Fig. [16] (upper 
panel) shows the initial setup for all the runs, which consist of a cir- 
cular cylinder plus a Cartesian background grid of 250 x 50 mesh 
generating points. The dynamic viscosity coefficient p takes five 
different values: 2.5 x 10" 2 , 5 x 10" 3 , 2.5 x 10" 3 , 1.25 x 10" 3 
and 8.3 x 10 -4 . These values correspond to Reynolds numbers of 
2, 10, 20, 40 and 60. 

For each one of the tests, we show the resulting streamlines at 
time t = 9.9 (or an equivalent dimensionless time of i = tU/d « 
31.0) in Fig. [17] Below Re ~ 40, the flow is steady and symmetric 
above and below the cylinder. As the Reynolds number increases, 
the size of the wake behind the cylinder grows. Although in this 
example the structure of the wake is poorly resolved, the increase 
in Re is accompanied by an increase of vorticity confined within 
the wake. 

Above Re ~ 40, the wake behind the cylinder starts to be- 
come unstable. This can be clearly seen in the streamline pattern of 
the Re = 60 panel. As the wake becomes unstable, the symmetry 
between the upper and lower portions of the domain is broken, at 
which point the flow becomes unsteady, such that the streamlines 
are no longer a valid representation of the Lagrangian trajectories 
of fluid parcels. This marks the onset of the von Karman vortex 
street, and the eventual transition to fully developed turbulence. 

To further illustrate the flexibility of the mesh construction in 
AREPO, we can repeat this experiment with the mesh generating 
points set to remain static, thus recovering an Eulerian grid code. 
In addition, we increase the resolution by a factor of four to bet- 
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Figure 15. Geometry of our wind tunnel set-up with a circular obstacle. 





Figure 16. Mesh near a circular cylinder inside a wind tunnel. The mesh 
contains both stationary mesh-generating points (defining the solid cylin- 
der and two layers of cells used to create the cylindrical solid surface) 
and moving mesh-generating points (the remainder of the grid). The up- 
per panel shows the initial setup, which highlights the cells representing 
the solid cylinder, and the two layers of fluid cells for which the equations 
of hydrodynamics are solved as in a standard stationary mesh. The total 
number of cells in the wind tunnel is 12,478 (roughly 250 x 50). The 
perimeter of the cylinder is outlined by 30 cells, and its diameter is equiv- 
alent to eight cells across. The Voronoi faces in between the red and grey 
cells define the boundary at which the no-slip condition is imposed. The 
lower panel shows the same region at a later time. Whereas one layer of 
mesh-generating points surrounding the cylinder has remained stationary, 
the rest of the background mesh has moved downstream and transformed to 
a generic unstructured Voronoi mesh as it moves along with the fluid. 



ter resolve the wake behind the cylinder. In Fig. [18] we show the 
density contrast for five different Reynolds numbers. For station- 
ary flow, the density distribution traces the streamline topology. At 
Re = 100, we can appreciate how the fully developed von Karman 
vortex street looks for a compressible gas. 



4.6 Three Dimensions: Taylor- Couette Flow 

Circular Couette flow is a stable, special case of the more complex 
and richer three-dimensional Taylor-Couette flow (Taylor 1923). 
Taylor found that when the angular velocity of the inner cylinder is 
increased above a certain threshold, Couette flow becomes unsta- 
ble. After this transition, different states have been identified, the 
most famous of which is the Taylor vortex flow, characterized by 
axisymmetric toroidal vortices. The diversity of states for Taylor- 
Couette flow has been explored in the past, most notably by Coles 
( |1965| ) and |Andereck et al.|(T986] ). The latter work lists up to 18 dif- 
ferent flow regimes observed in flow between independently rotat- 
ing cylinders. Its ' Andereck diagram", which explores the stability 
of the Taylor-Couette problem for a variety of Reynolds numbers, 
has become the standard benchmark for computational experiments 
of flow between rotating cylinders. 

Although the computational and experimental study of three- 
dimensional Couette flow peaked during the 1980' s with the classi- 
cal works of Anderec k et al.| ( p"986| ) andJMarcus (1984a. b), in recent 
years it has regained popularity (e.g. Dong 2007 ; Avila et al. 2008 , 
Meseguer et al. 2009a b ) mainly driven by the experimental studies 
of magnetized and unmagnetized rotating flows of |Ji et al.| (2001 
2006) and Sisan et al. (2004), which have resulted in significant 
progress on the characterization of the magnetorotaional instability 
(MRI; Ba lbus & Hawley|1998| ) in the laboratory. 

In this section, we briefly explore the evolution of Taylor- 
Couette flow on a moving Voronoi mesh. Although the AREPO code 
is not specifically designed for problems with symmetric geome- 
tries where static cylindrical meshes have proven to be more suit- 
able, we have included this test to emphasize that our method works 
in three dimensions in an analogous way to the two-dimensional 
examples shown above. It is straightforward to extend the two- 
dimensional Couette flow shown above to three dimensions using 
AREPO. Since the mesh is obtained from a distribution of mesh- 
generating points, all that is needed is to replicate the initial condi- 
tions shown in Figure [13] in the vertical direction (about 80 times) 
to fill up a cubic box. 
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(a) Re = 2 




(b) Re = 10 




(c) Re = 20 




(d) Re = 40 




(e) Re = 60 
Figure 17. Streamlines for compressible gas flow around a cylinder at five different Reynolds numbers, as labeled. 



A standard validation for a Taylor-Couette simulation with 
azimuthal and axial periodicity may include, for example (Avila, 
private communication): obtaining perfect axial symmetry at low 
Reynolds number (i.e. circular Couette flow ), followed by obtain- 
ing the first bifurcation to axially symmetric Taylor vortices, and 
by reaching the second bifurcation to wavy vortices. These transi- 
tions occur sequentially as the angular velocity of the inner cylin- 
der is increased while keeping the outer cylinder stationary (see the 
phase diagram of Andereck et al. 1986 ). However, it is not the pur- 
pose of this section to explore these transitions exhaustively; we 
only want to show that the third dimension works with our tech- 



nique. We thus have focused on a particular configuration: counter- 
rotating Taylor-Couette flow, for which it is easy to obtain axially 
symmetric Taylor vortices (although these might relax back to Cou- 
ette flow after several rotation periods; e.g. |Liao et al.|19 99). The 
geometry described in Fig.[TT]is replicated in the vertical direction 
such that the computational domain is now a cube of dimensions 
6x6x6, with periodic boundary conditions in the z-direction. 
The initially Cartesian mesh will eventually relax in all directions 
as the flow evolves (Fig. [T9j. The cylinder is effectively infinite, 
like in the two-dimensional case, except that this time there is no 
imposed symmetry along the ^-direction. We choose the cylinder 
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Figure 18. Density contrast of compressible flow past a cylinder at five different Reynolds numbers, corresponding to Re = 2, 10, 20, 40, 100, from top to 
bottom. All five numerical experiments were computed with a static Cartesian mesh at moderately high resolution (1000 x 250), where the cell size is 1/32 of 
the cylinder's diameter. 
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Figure 19. Vertical slice of the three-dimensional Voronoi tessellation in 
Taylor-Couette flow at the time Taylor vortices have developed. This same 
slice is used when visualizing the v x ,v y and v z fields (Fig. [20]) and Fig.|2T). 



radii to be Ri = 1.0 and R2 — 2.5, just like in the 2D exam- 
ple, and the respective angular velocities are Q± = 0.8 ^2 = —0.5 
(counterrotating). The dynamic viscosity is \i — 0.005 and the fluid 
is started from rest with p — P — 1. The inner and outer Reynolds 
numbers (Rei = i?A(i?2 — Ri)p//jl; e.g. Liao et al. 1999) are 
Ri = 240 and R2 = —375, respectively. 



The geometry of the problem is shown in Fig. [20ft. A vertical 
slice is taken at a time when the Taylor vortices have developed 
(the corresponding sliced mesh is shown in Fig.[l9]). The azimuthal 
velocity on that slice shows deviations from the symmetry in z 
present in the circular Couette regime (Fig. |20]p). Looking at the 
poloidal velocity field on that same slice (v x and v z in Fig.|2TJ one 
can appreciate, near the inner cylinder, the circular vertical motion 
characteristic of the Taylor vortices. 

In Figure|2T] we show the velocity field of this Taylor-Couette 
experiment at time t = 128 for two different slices of the vol- 
ume: (a) along the x-axis, and (b) along the y-axis (i.e. at 90° from 
the first slice). Except for the numerical noise, the two solutions 
are nearly indistinguishable, evidence of a global axially symmet- 
ric Taylor vortex flow (for a very similar configuration, see Fig. 3 
in |Liao et al.|1999| >. This flow starts to develop at time t ~ 60 and 
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remains essentially unaltered for several rotation periods. At much 
longer time scales, the flow would presumably decay back to a two- 
dimensional Couette flow as seen in the roughly similar test carried 
out by |Liao et al.]fT999] ). 



5 CONCLUDING REMARKS 

We have presented a new numerical approach for solving the three- 
dimensional, compressible NS equations on a dynamic mesh using 
the new astrophysical code AREPO. This novel approach, an exten- 
sion of the finite volume method, defines the computational mesh 
as a Voronoi tessellation moving with the local flow. The advan- 
tages of using a dynamic Voronoi mesh for transient and station- 
ary flows under diverse boundary conditions has been addressed. 
The implicit adaptivity of the quasi-Lagrangian mesh elements, in 
addition to the well-behaved topological properties of the Voronoi 
tessellation, ensure both geometric flexibility and low numerical 
diffusivity. In addition, the shock capturing, second-order-accurate 
finite-volume scheme implemented in the rest-frame of each mov- 
ing cell provides high accuracy. 

We have described in detail the algorithm used to estimate the 
viscous diffusion of momentum across inter-cell boundaries. Our 
scheme produces smoothly varying estimates of the viscous terms, 
resulting in accurate and stable solutions. The method extends pre- 
viously known finite-volume formulations of the NS equations with 
the introduction of a new reconstruction scheme that represents a 
compromise between the use of piece- wise constant gradients and 
fully consistent quadratic-reconstruction schemes. 

For pure hydrodynamic flow, the CPU time consumption of 
our code per timestep is typically quite a bit higher than for struc- 
tured mesh codes or SPH codes, for the same number of resolution 
elements. In three dimensions, the factor is close to 2 relative to 
SPH (if 64 smoothing neighbors are used), and up to ~ 3 rela- 
tive to a Cartesian mesh codes. The additional computational time 
goes mostly into the Voronoi mesh construction overhead, which 
is simply not needed by a structured mesh code, and also into an 
enlarged computational cost for the flux computations. The latter 
comes about because of a larger number of faces per cell (in 3D, 
there are 6 sides for a cubical cell, but for a 3D Voronoi mesh, 
we have of order ^12 sides for each polyhedral cell). It is how- 
ever important to note that other, problem-dependent factors should 
be taken into account when assessing the performance in practice. 
For example, if large bulk velocities are present, our method can 
take considerably larger timesteps than a corresponding fixed mesh 
code. Also, because the advection errors are reduced in our scheme, 
fewer cells are required to reach a given accuracy, so that our code 
can then end up being computationally more efficient. We also note 
that once self-gravity is included (as in many of our primary tar- 
get applications in astrophysics), the relative speed difference in 
the hydrodynamic part between the structured fixed mesh and our 
moving Voronoi mesh becomes much less of an issue, because the 
cost of calculating self-gravity sufficiently accurately for arbitrary 
geometries substantially reduces the relative importance of the hy- 
drodynamical cost. 

As part of our study, we have verified the reliability of our 
new method through a series of example calculations that range 
from simple flows with known analytic solutions to traditional ex- 
periments of well-known quantitative behavior. The demonstrated 
ability of the scheme to reproduce exact solutions as a function of 
time, even if the flow is started impulsively from rest, is reassur- 
ing. These examples also show the flexibility of the scheme in the 



presence of different solid surfaces moving in diverse ways. In all 
of these examples, the overall structure of the numerical scheme 
is identical, and the boundary conditions are set solely by the pre- 
scribed motion of the surfaces, which consist of controlled collec- 
tions of Voronoi cells. 

Although we have tested the performance of AREPO in con- 
figurations possessing a high degree of symmetry, it is in complex 
asymmetric problems where the moving-mesh approach would 
show all its power. The flexibility of the Lagrangian nature of 
the mesh will allow us to simulate, for example, complex astro- 
physical objects where viscosity is presumed to play a signifi- 
cant role. One such problem is the simulation of accretion disks 
around young stars. Although angular momentum transport in ac- 
cretion disks is attributed to turbulence (most likely of magneto- 
hydrodynamic nature), this process is usually modeled both an- 
alytically (e.g. Shakura & Sunyaev 1973, Lynden-Bell & Pringle 
^974l|Pringle|1981HLin & Pringle|1987) as well numerically (e.g. 



Kley & Lin 1992 ; Masset 2000 ; D' Angelo et al. 2002 de Val-Borro 



et al. 2006, Paardekooper & Mellema 2006, Mudryk & Murray 



2009 just to name a few) by laminar flow in the presence of turbu- 
lent viscosity (Boussinesq approximation to eddy viscosity), due to 
the computational cost of global models of magneto-hydrodynamic 
disks. Another application of viscous flow is the plasma viscosity 
at galaxy cluster scales (e.g. Sijacki & Springel 2006). However, it 
is likely that in such systems viscosity, as well as thermal conduc- 
tion, is anisotropic (Braginskii 1965] see |Dong"& Stone 2009] for 
an example). In such a case, the viscous stress tensor in Eq. ([6} can 
be easily generalized to include the up to seven independent vis- 
cosity coefficients (Lifshitz & Pitaevskii 1981). It will be particu- 
larly exciting to couple the local anisotropy directly to the magnetic 
field topology, with the latter calculated self-consistently using a 
recent magnetohydrodynamics implemention in AREPO (Pakmor 
|etal,|2011| . 

Its powerful flexibility will make AREPO an interesting code 
both for astrophysical simulations of viscous flow, but potentially 
also in engineering applications where the ability to cope with 
curved and moving boundaries is particularly attractive. 
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v y at 2=3 





(a) 



(b) 



Figure 20. a) Illustration of the three-dimensional flow between two independently rotating cylinders. The figure shows a plane along the radial direction 
where the local velocity field is evaluated, b) Velocity field in the y-direction between the cylinders for a slice defined by y = 3 (i.e. along the diameter 
of both cylinders). For this particular plane, v y is equivalent to the azimuthal velocity vq. The color scale goes from v y = Q2R2 = —1.25 (blue) to 
v y = Q±Ri = 0.8 (red). This example shows that vq is no longer independent of z. Thus the two-dimensional solution of Eq. (65) is no longer valid. 
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Figure 21. Velocity structure of axisymmetric Taylor vortex flow at time t = 119 at two different meridian planes separated by 90°. The poloidal velocity 
field (vr, Vz) is color mapped in the linear range [—0.08 (blue), +0.08 (red)], while the azimuthal field (vq) is color mapped in the linear range [—1.2 (blue), 
0.6 (red)]. The streamlines illustrate the vector field in the poloidal plane, showing with clarity the nature of Taylor vortices. 
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APPENDIX A: GRADIENT EXTRAPOLATION 
COEFFICIENTS 

The extrapolation of the velocity gradients (e.g. Eq.[36} requires a 
numerical estimate of the gradient matrix as well as an estimate for 
the time derivative of the gradient. For the latter, the tensors A a pb 
and B a pba are needed (Eq. [34l. Both tensors depend on the cell- 
centered scalar quantities as well as their gradients. The values of 
A a f3b are (e.g. Toro 2009) 
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The tensor B ai3ba = d a A a(3b = A aj3b , a (with a, b 
x, y, z or 1, 2, 3 and a, f3 = 0, 1, 2, 3, 4) has components: 
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